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I.  Introduction 


We  have  successfully  completed  our  first  version  of  an  ab-initio  molecular 
dynamics  computer  code<^whicfi*  can  simulate  the  motion  of  atoms  at  a  surface  and 


in  the  bulk  of  a  semiconductor.  We  combineJNewton's  equation  for  the  nuclei, 

i  r  r  p  ’■ '  ‘ 

F  =  ma^v  with  the  Schrodinger  equation  for  the  electrons,  HT=E $  to  obtain  a 
uniform  picture  of  a  covalent  systems  dynamical  properties. 

,,  per  ruv-AW 

'-Tn  this  document,  we' will  list  some  of  the  simulations , that  we  have 

L  'j 

performed.  These  simulations  should  not  be  considered  as^being  final,  but  we 
decided  to  try  a  variety  of  problems  without  going  into  extreme  depth,  so  that 

we  can  better  assess  the  strengths  and  weaknesses  of  the  technique. 

/  .v  ji>  "i- 

We  have  developed^ a  tight-binding  method^ -whtrse^tight^binding  matrix 
elements  acvt' calculated  entirely  from  first  principles.  No  fitting  to 
experiment  of  any  quantities  is  needed  or  done.  We-u se this  tight-binding 
Hamiltonian  calculate  the  electronic  structure  of  the  material.  The 
electronic  structure  theory  is  based  on  the  local  density  approximation  (LDA), 
with  no  adjustable  parameters .~^We  have  made  a  number  of  approximations  which 


are  detailed  in  Refs.  1,  2,  and  3.  The  use  of  approximations  is  essential  to 
having  a  method  which  is  fast  enough  to  be  useful  for  simulation  of  medium  and 
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large  size  systems.  Where  possible,  ve-have-rompared  oujt  results  to 
experimental  data  and  find  agreement  consistent  with  the  LDA.  '  •  ... 
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We  calculate  forces  on  the  Si  atoms  by  employing ' the  Hellmann-Feynman 


theorem.  As  a  consequence  of  our  ab-initio  tight  binding  approach,  in 
conjunction  with  a  careful  optimization  of  computational  techniques,  ve  find 
that  it  is  practically  possible  to  perform  molecular  dynamics  simulations  in 
supercells  containing  up  to  a  couple  hundred  atoms.  We  are  working  on  dramatic 
generalizations  and  enhancements  of  the  procedure,  and  expect  to  make  rapid 
progress. 

We  close  this  section  with  a  mention  of  some  general  features  of  the 
simulations  we  have  performed.  In  each  simulation  we  choose  initial  conditions 
so  that  the  linear  and  angular  momentum  is  zero,  and  the  center  of  mass  is 
stationary.  Additionally,  we  will  always  use  the  convention  that  "temperature" 
T  will  refer  to  "kinetic  temperature",  that  the  value  of  T  such  that  (3/2)kT  = 
average  kinetic  energy  per  particle  of  the  system. 

The  rest  of  this  report  gives  brief  examples  of  a  variety  of  simulations. 


II.  Brief  Examples  of  Simulations 

For  some  simulations  we  have  produced  a  plot  of  "kinetic  temperature"  T  vs. 
time,  where  T  satisfies  the  equation:  average  kinetic  energy  per  particle  = 
(3/2)kjjT.  This  type  of  plot  illustrates  the  time  dependence  of  the  kinetic 
energy  of  the  system.  As  we  expect  from  classical  equipartition,  a  system 
which  starts  from  equilibrium  at  temperature  2T  rapidly  converts  about  half  of 
that  energy  into  potential  energy.  The  temperature  then  fluctuates  about  T. 
Thus,  for  a  system  started  with  equilibrium  geometry,  we  may  confidently 
predict  that  the  final  average  temperature  will  be  quite  close  to  half  of  the 
initial  kinetic  temperature  we  started  with.  Note  that  this  comment  is  valid 
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only  for  cases  where  we  start  from  eqilibrium,  however. 
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A  second  kind  of  plot  we  provide  is  what  we  term  a  "snapshot  graph",  which 
for  planar  motions  gives  the  trajectories  of  each  particle  in  the  central 
supercell.  For  non-planar  dynamics,  we  show  a  projection  onto  the  xy  plane. 

Another  plot  which  lends  insight  into  the  molecular  dynamics  is  a  graph  of 
<r^(t)>,  which  illustrates  how  far  an  "average"  atom  deviates  from  its  starting 
position  at  time  t  from  where  it  began  at  t=0.  Here  <>  is  a  configuration 
average  over  all  particles  in  a  cell. 

We  have  in  many  cases  performed  a  "dynamical  quenching"  of  the  system.  By 
this  we  refer  to  a  procedure  designed  to  find  a  geometrical  configuration  with 
minimum  potential  energy.  The  minimum  found  may  be  a  global  minimum  (as  we 
will  illustrate  for  small  clusters  of  atoms),  or  a  local  minimum  for  a  more 
complicated  example,  such  as  the  dynamical  quenching  of  a  large  number  of  atoms 
at  a  very  high  temperature.  The  method  we  found  most  effective  was  to  quench 
the  atomic  velocities  (i.e.  set  them  to  zero)  near  a  peak  of  the  temperature 
vs.  time  plot,  then  let  the  system  evolve  naturally  until  it  reached  the  next 
maximum,  quench  again,  and  so  on.  This  was  found  to  be  a  very  effective  and 
accurate  means  of  finding  minimum  energy  configurations. 

For  some  of  our  calculations  we  found  it  useful  to  calculate  double  time 
autocorrelation  functions  of  two  kinds.  These  are  convenient  because  they 
probe  the  "memory"  of  the  system.  Roughly  speaking,  we  may  interpret  small  (in 
some  appropriate  sense)  values  of  an  autocorrelation  function  evaluated  at  two 
different  times  to  mean  that  the  dynamical  variable  from  which  we  construct  the 
correlation  function  is  only  weakly  correlated,  or  entirely  uncorrelated.  We 
may  think  of  this  as  a  gauge  of  "information  loss"  with  time  evolution.  Some 
specific  functions  we  calculate: 

1.  A  velocity- velocity  autocorrelation  function: 

«(t)  -  tfT)-:(0)>  , 

<v(0)-v(0)> 


where  a  cell  average  <>  is  taken  of  the  dot  product  df  atomic  velocities  at  two 
times  t}  and  t2»  and  where  x  =  t2  -  tj.  Since  this  function  is  invariant  under 
translations  of  the  temporal  origin,  we  may  take  g  to  be  a  function  only  of  T. 
2.  A  position-position  autocorrelation  function: 

Gx(t)  =  <u( t) • u(0) >  , 

where  x  is  the  same  as  for  the  previous  equation.  Here  u  is  the  departure  from 
equilibrium.  Function  gx  obviously  contains  information  about  the  "diffusive" 
properties  of  the  molecular  motion. 

In  either  case,  it  is  worthwhile  to  further  consider  the  Fourier  transforms 
of  the  correlation  functions.  The  frequency  domain  version  of  g  and  gx  is 
clearly  the  natural  function  to  work  with  if  we  are  interested  in  a  normal  mode 
view  of  the  motion.  In  this  case,  the  frequency  functions  are  interpreted  as  a 
spectral  mode-density. 

II-A  Si  Clusters. 

The  simplest  class  of  problems  susceptible  to  our  techniques  is  the 
dynamics  of  small  silicon  clusters.  Our  technique  does  not  require 
periodicity,  so  that  reciprocal  lattice  vectors  G  are  never  used.  This  allows 
molecules,  bulk  systems,  and  surfaces  to  be  handled  all  within  the  same 
framework. 

For  Si  clusters  we  investigate  vibrational  modes,  equilibrium  ground  state 
geometries,  electronic  states,  high  temperature  phenomena,  and  collisions.  The 
following  simulations  illustrate  these. 

Six  at  room  temperature 

An  Sij  molecule  is  started  at  equilibrium  with  random  planar  velocities 
having  an  average  energy  of  600K.  The  center  of  mass  is  at  rest  and  the 
molecule  has  no  angular  momentum.  This  constrains  the  motion  to  two 


dimensions. 


A  snapshot  of  the  simulation  shoving  the  positions  of  the  three  atoms  in 
the  Si3  cluster  over  the  time  of  the  simulation.  The  motion  is  entirely  in  the 
plane.  There  are  2200  time  steps  each  of  1.57  fs. 

Figure  2. 

The  kinetic  "temperature"  as  a  function  of  time  for  the  Si3  simulation. 

The  system  starts  at  600K,  but  the  average  over  the  whole  simulation  is  close 
to  300K  as  expected  from  the  equipartition  theorem.  Notice  that  the 
"temperature"  of  the  molecule  oscillates  between  600K  and  OK. 

Si^  at  -2500K 

We  now  consider  what  happens  to  an  Si3  molecule  at  high  temperature.  A  Si3 
molecule  is  started  with  random  planar  velocities  having  an  initial  average 
energy  of  3000K.  The  atoms  were  not  started  from  equilibrium;  rather  they  were 
started  from  a  right  isosceles  triangle  with  a  bond  length  of  2.0  A.  The 
average  temperature  over  the  entire  simulation  is  -2500K.  The  center  of  mass 
is  at  rest  and  the  molecule  has  no  angular  momentum.  This  constrains  the 
motion  to  two  dimensions.  There  were  3600  time  steps,  with  each  time  step 
being  1.57  fs. 

Figure  3. 

A  snapshot  of  the  time  evolution  of  the  Si3  simulation  at  -2500K.  The 
motion  is  entirely  2  dimensional,  and  is  very  far  from  harmonic. 

Figure  4. 

The  temperature  of  the  system  as  a  function  of  time.  The  system  started 


with  an  average  kinetic  energy  of  3000K,  but  higher  temperatures  are  seen  since 
the  system  was  not  started  from  the  equilibrium  geometry. 


A  dynamical  quenching  was  performed  on  an  Sij  molecule  to  find  the  ground 
state  structure.  The  Si3  molecule  was  started  in  a  right  triangle  with  sides 
of  lA  and  2k.  The  atoms  start  very  close  together,  and  in  fact  the  atoms  would 
fly  apart  if  they  were  not  annealed.  However,  since  it  is  quenched,  the  system 
relaxes  rapidly  to  the  equilibrium  structure  in  -100  time  steps.  The  final 
geometry  is  thk.  of  an  isosceles  triangle. 

Figure  5. 

A  snapshot  of  the  dynamical  quenching  of  the  Si3  molecule.  The  final 
positions  are  indicated  by  the  ''dots”. 

Figure  6. 

The  temperature  as  a  function  of  time  during  the  quench.  We  use  the 
annealing  procedure  discussed  in  Section  II. 

Figure  7. 

The  final  minimum  energy  geometry  of  the  Si3  molecule. 

Figure  8. 

The  energy  eigenvalues  of  the  single  particle  LDA  Hamiltonian,  as  a 
function  of  angle  for  the  Si3  molecule.  The  bond  length  is  kept  fixed  at 
2.189A.  The  Fermi  energy  shown  in  the  figure  separates  occupied  from 
unoccupied  states.  For  the  case  of  an  equilateral  triangle  (60  degrees),  the 
system  has  partially  occupied  levels,  and  is  unstable. 

Vibrational  spectra  of  Si-^ 

Figure  9. 

The  spectrum  of  vibrational  modes  obtained  from  the  «-transform  of  the 


velocity-velocity  autocorrelation  function,  <  v(x)v(0)  >/<  v(0)v(0)  >.  The 
peaks  occurred  at  normal  mode  frequencies,  with  broadening  determined  by  the 
total  elapsed  time  of  the  simulation.  (In  this  figure,  3600  steps  of  1.57  fs 
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were  used).  Two  different  average  temperatures  are  shown.  The  simulation  has 
no  angular  momentum  or  center  of  mass  motion  leaving  the  number  of  vibrational 
modes  at  3n-6. 

Figure  10. 

We  show  for  completeness  the  equilibrium  structures  determined  through 
dynamical  quenching  for  Si4,  Si5,  Sig,  and  Siy. 

Collisions  involving  five  Si  atoms. 

Our  final  simulation  involving  Si  clusters  involves  the  collision  of  two 
molecules  to  form  an  Si^  molecule. 


The  simulation  begins  with  two  Si2  molecules  coming  together  and  colliding 
to  form  Si4*  The  center  of  mass  motion  and  total  angular  momentum  are  zero. 
When  the  two  collide,  a  large  energy  of  bonding  is  converted  into  kinetic 
energy  and  the  newly  created  Si4  molecule  vibrates  violently.  The  vibration 
continues  for  a  time,  and  finally  a  dynamical  is  performed  to  bring  the  system 
to  its  equilibrium  geometry. 

Figure  11. 

The  snapshot  of  the  time  evolution  of  the  collision  between  two  Si2 
molecules  to  form  a  highly  excited  Si4  molecule. 

Figure  12. 

The  temperature  as  a  function  of  time  for  the  collision  of  two  Si2 
molecules.  The  first  1000  time  steps  are  free  motion  of  non-interacting  Si2 
molecules.  At  time  step  -1000  the  two  molecules  interact  and  the  kinetic 
temperature  rises  substantially.  The  newly  formed  Si4  molecule  evolves 
naturally  in  time  until  about  time  step  2500  when  the  system  is  dynamical 
quenched  equilibrium. 
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[100]  phonons  in  bulk  Si. 

We  now  compute  dynamically  the  phonon  vibrational  modes  of  bulk  Si  for 
phonons  in  the  [100]  direction.  A  supercell  is  used  which  contains  16  Si 
atoms.  Eight  special  k  points  are  used  from  the  entire  Brillouin  zone  in  the 
computation  of  the  total  energy  from  which  the  forces  are  determined. 

Figure  13. 

(a)  The  [0  -1  1]  plane  in  the  diamond  lattice,  from  which  the  supercell  is 
constructed. 

(b)  The  supercell  geometry  used  to  determine  phonon  modes  with  wavevector 
along  X.  The  cell  contains  four  segments  of  the  type  shown  in  (a).  The  larger 
atoms  are  shifted  above  the  plane  by  a//2,  where  a  is  the  cubic  lattice 
constant  of  the  diamond  lattice.  There  are  16  atoms  in  this  cell  (some  edge 
atoms  are  in  different  cells),  and  the  atoms  in  this  cell  can  vibrate  with 
phonon  modes  of  wavevector  0/4,  1/4,  2/4,  3/4,  and  4/4  of  the  X  point  in  the 
Brillouin  zone. 

Simulation  of  all  modes  with  k  along  [100]. 

This  simulation  shows  bulk  phonon  modes  of  all  polarizations  with 
wavevector  along  [100].  The  system  is  started  from  the  equilibrium 
configuration,  and  the  atoms  are  given  random  velocities  with  an  initial 
temperature  of  600K.  This  yields  an  average  over  the  entire  simulation  of 
-300K.  The  motion  consists  of  a  linear  combination  phonon  of  the  five 

Vavevectors  (0,  0.2,  0.4,  0.6,  0.8,  and  1.0  of  —  (100)),  with  all 

Si 

polarizations. 

Figure  14. 

A  snapshot  of  the  first  few  atoms  in  the  supercell  for  the  simulation  with 
all  modes  present.  Notice  that  for  k  along  [100],  the  amplitude  is  much 


greater  perpendicular  to  k,  than  parallel  to  k.  This  is  due  to  the  low- 
frequency  transverse  acoustic  mode. 

Figure  15. 

The  double  fourier  transform  of  the  transverse  velocity-velocity 
autocorrelation  function,  gn(t),  where  n  is  an  index  for  particle  n.  The 
transform  to  co-space  was  described  in  the  introduction.  The  new  feature  of  the 
double  transform  is  the  use  of  a  spatial  k  transform  with  the  atomic  position 
of  atom  n.  Thus  gn(t)  is  transformed  to  g(k,co).  This  function  is  plotted  as  a 
function  of  «  for  the  various  values  of  k.  Turning  the  figure  on  its  side 
yields  the  more  conventional  w  vs.  k  dispersion  relation. 

Figure  16. 

Conventional  phonon  dispersion  relation  calculated  by  computing  the 
dynamical  matrix  using  the  frozen  phonon  technique  in  the  supercell  geometry 
describe  above.  The  atoms  are  moved  by  small  amounts  statistically  from 
equilibrium,  and  effective  spring  constants  are  determined.  About  100  k 
vectors  in  the  Brillouin  zone  are  used  to  evaluate  the  electronic  total  energy. 
The  theory  is  in  close  agreement  with  experiment. 

Simulation  of  lattice  relatation  around  Si  vacancy. 

We  have  simulated  the  lattice  relaxation  around  a  Si  vacancy.  We  use  a  32 
atom  BCC  supercell  geometry  for  the  perfect  crystal.  Since  one  site  is  a 
vacancy,  the  atom  at  the  origin  is  missing  and  there  are  actually  only  31 
atoms.  The  lattice  constant  chosen  is  5.5A,  as  this  corresponds  to  the 
theoretical  bulk  minimum.  Four  special  k-points  are  used.  They  are 
0.285593(1,1,1),  0.285593(-l,-l, 1),  0. 285593(-l, 1,-1) ,  and  0.285593(1,-1,-1)  in 
inverse  A  units.  These  are  the  Monkhorst-Pack  special  points  for  the  BCC 
lattice  with  lattice  constant  5.5A. 
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The  simulation  starts  at  the  "ideal"  positions,  meaning  the  atoms  are  in 
the  same  position  as  if  the  vacancy  were  not  there.  We  give  the  atoms  random 
velocities  with  average  temperature  of  600K.  Normally  this  would  yield  after 
some  time  an  average  velocity  of  300K  (room  temperature).  After  200  time  steps 
(each  of  .206  fs),  we  begin  to  quench  the  system. 

Figure  17. 

The  electronic  energy  level  spectrum  in  the  band  gap  region  for  the 
vacancy.  A  triply  degenerate  level  is  found  in  the  band  gap  which  is  occupied 
by  two  electrons.  This  corresponds  to  the  situation  at  time  t*0  (unrelaxed). 
Because  the  levels  are  degenerate  and  partially  occupied,  the  system  can 
undergo  a  Jahn-Teller  distortion  in  which  the  symmetry  is  broken.  The  level 
structure  at  the  end  of  the  anneal  is  shown  on  the  right.  A  single  level  is 
found  at  lower  energy  and  is  doubly  occupied.  A  pair  of  doubly  degenerate 
levels  are  found  at  higher  energy  are  are  unoccupied  (for  the  case  of  the 
neutral  vacancy). 

Figure  18. 

The  final  positions  of  the  first  nearest  neighbors  around  the  vacancy  drawn 
to  scale.  The  square  shows  the  geometry  of  the  "ideal"  unrelaxed  system,  and 
the  rectangle  the  geometry  of  the  Jahn-Teller  distorted  system.  The  Jahn- 
Teller  distortion  is  of  tetragonal  symmetry  with  the  x  and  z  displacements 
being  equal,  while  the  y  displacement  is  different.  The  final  tetragonal 
pattern  could  have  been  either  x-,  y-,  or  z-like.  It  came  out  y-like  in  this 
simulation,  because  there  was  no  energy  barrier  from  the  configuration  which 
existed  when  the  annealing  began  at  the  200th  step. 

Free  vibrations  and  simulated  annealing  of  "liquid"  Si  at  1500K. 


In  this  simulation,  we  start  the  particles  at  equilibrium,  but  with  the 
very  high  average  velocity  of  3000K.  This  ought  to  yield  (according  to  the 


equipartition  theorem)  an  average  temperature  near  equilibrium  of  about  1500K. 
This  is  near  the  melting  point  of  bulk  Si  (expt.  1683K).  This  high  temperature 
stipulation  was  exploratory,  so  that  only  the  k  equal  to  zero  vavevector  is 
used  as  a  "special"  point  to  sum  over  the  Brillouin  zone.  Four  special  points 
should  be  used  for  a  more  rigorous  simulation. 

The  simulation  is  in  two  parts.  In  the  first  500  times  steps  (2.09  fs 
each),  the  particles  are  allowed  to  move  freely  under  the  equations  of  motion. 
They  can  vibrate  and/or  diffuse.  In  the  next  800  time  steps  the  system  is 
annealed  by  quenching  the  system  by  setting  the  velocities  equal  to  zero  on 
occasion. 


Figure  19. 

The  snapshot  of  the  first  500  time  steps  of  the  simulation  where  the  atoms 
are  moving  freely.  The  geometry  is  the  32  atoms  BCC  geometry.  One  notices 
large  excursions  from  the  equilibrium  positions,  indicating  diffusive  motion 
and  or  melting. 

Figure  20. 

The  average  "temperature"  as  a  function  of  time.  At  the  first  time  step, 
the  particles  are  started  from  equilibrium  with  an  average  velocity 
corresponding  to  a  temperature  of  3000K.  Within  a  very  few  time  steps,  the 
kinetic  energy  is  reduced  by  about  a  factor  of  2,  having  been  converted  into 
potential  energy.  The  average  temperature  over  the  whole  simulation  is  1400K, 
but  large  fluctuations  are  evident. 


Summary 

Ve  have  just  illustrated  a  few  applications  of  our  newly  developed  quantum 
molecular  dynamics  method  and  computer  code.  The  simulations  reported  here  are 
our  first  and  ve  anticipate  large  leaps  are  still  possible.  Ve  believe  our 


work  is  a  major  breakthrough  in  the  area  of  quantum  molecular  dynamics. 
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